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ABSTRACT 
Wetland's values are considered an important part of global ecosystems. The primary aim of this study was to develop a 
methodology for the assessment of wetland suitability potential by employing a Geographic information system (GIS) and Remote 
sensing (RS) technique. This method applied to the Federal University of Agriculture Abeokuta, Ogun State, Nigeria. Soil samples (0 
to 30cm depth) were collected, analyzed, and interpolated using Inverse distance weighted techniques in ArcGIS 10. Supervised 
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classification algorithm was used to get land cover information using Landsat 8 Operational land images. A GlS-based land 


suitability analysis was performed after identifying the following criteria: Land surface temperature, Modified normalized difference 
water index, Land use/land cover, Soil adjusted vegetation index, and Enhanced built-up and bare land index of the study area. The 
relative importance of each criterion was assessed utilizing the Analytical hierarchy process in Multi-criteria decision making. Land 
suitability index was calculated and a suitability map was produced. The result of the suitability map shows that 28.9% and 33.2% 
was moderately and marginally suitable, while 37.9% was marginally not suitable for rice production. The present study 
demonstrated the efficiency of GIS and RS data for wetland sustainability. 


Keywords: Wetland; rice; geographic information system; remote sensing; analytical hierarchy process. 


1. INTRODUCTION 


Wetlands are areas whose formation is influenced by ecological, hydrological, and geomorphological processes (Frenken and 
Mharapara, 2002). The desire to find the right balance between the exploitation of wetland services and their conservation has led to 
an increase in wetland benefits to humans (Mitsch and Gosselink 1993). Wetland degradation may increase in the next 50 years if 
wetlands ecosystem services and functions are not effectively managed and monitored (MA, 2005b). Moreover, due to the 
increasing world population and the growing pressure on food production, sub-Saharan Africa such as Nigeriamay experience food 
scarcity by the year 2020 (Joshua et al. 2013). In recent decades, the use of wetlands has increased significantly in many developing 
countries (McCartney et al. 2010). In Nigeria, wetland areas cover 7.2 percent of the total land area of the country (Ojanuga, 2006). 
However, the occurrence of wetland soils in Nigeria has been associated more with rice production with three landforms namely; 
inland depressions, floodplains, and coastal plains (Fasina, 2005). 

Rice (Oryzasativa) belongs to the Gramineae family and constitutes a significant component of major staple food crops 
(Ayolagha et al. 2012). Suitability of wetland soils for optimum rice production depends on the underlying parent material, 
topography, land use, temperature, and rainfalletc. (Ahukaemere and Akpan, 2012). However, many wetland soils around the world 
are becoming more unsustainable. Hence, calls for geospatial techniques that serve as a bottom-up approach for consistent 
mapping at national, regional, and local levels (Karimi and Zeinivand, 2019). Land suitability analysis is a geographic information 
system (GIS) based application used in sustainably determining crop suitability (Zubair, 2006; Mustafa, 2011; Saaty, 1977; Jafari and 
Zaredar 2010; Vargahan et al. 2011). On the other hand, Remote Sensing (RS) techniques serve as valuable tools for the assessment 
of the well-being of wetlands at different scales (Merollaet al.1994; Davis, 2009). In the last few decades, integration of RS and GIS 
using the Analytical hierarchy process (AHP) has been widely applied to land suitability assessment problems (Kahinda et al. 2008; 
Mahmoud and Alazba 2014; Thapa and Murayama 2008; Cengiz and Akbulak 2009; Patil et al. 2012. The AHP serves as a powerful 
pair-wise comparison matrix that supports spatial decision system on the judgments of expert's priority scales (Thapa and 
Murayama 2007). 

In the present study, remote sensing data and a GlIS-based conceptual framework is applied with the Multi-criteria decision- 
making (MCDM) technique using AHP. The conceptual framework will help to identify the potential suitability of the wetland soils for 
rice production at the Federal University of Agriculture Abeokuta, Ogun State, Nigeria using geospatial techniques. 


Study Area 

The methodology of the land suitability evaluation was applied to the Federal University of Agriculture Abeokuta (FUNAAB) located 
next to the Ogun-Osun River Basin Development Authority (OORBDA) in the Northern part of Abeokuta city, Ogun state, Nigeria. 
The area is geographically described by latitudes 7° 13’N and 7° 20'N and by longitudes 3° 20’E and 3° 28’E enclosing approximately 
10,000 hectares (ha) of land area (Figure 1). It is characterized by undulating with extensively mild slopes bounded into six zones and 
punctuated in parts by ridges, isolated residual hills, plateaus, and valley landscapes with low lands (Savage, 2010). The area is 
mainly drained by the Ogun River and other streams namely: Oshinko, Ole, Alakata, Arakanga, Pala, Olu, Tigba, and Ajigbayin 
(Ufoegbune et al. 2010). 


Geology and soil 

The geology of the study area overlies metamorphic rocks of the basement complex, the majority of which are ancient being of pre- 
Cambrian age with great variation in grain size and mineral composition. This rock ranges from very coarse grain pegmatite to fine- 
grained schist and from acid quartzite to basic rocks consisting largely of amphibolites (Smyth and Montgomery, 1962). The area is 
sub-divided into sedentary, hill-creep, and hill-wash soils with the development of lowland forests (FAO/IUSS 2006; Ekanade, 2007). 
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The soils of the study area are dry for as long as 90 days and thus considered as the Udic moisture regime and therefore classified as 


Isohyperthermic temperature (Soil Survey Staff, 2010). For easy and efficient mapping of the study area, the soil classification map 
(Figure 2) produced by Sotona et al. (2013) was further modified and used for this study. 
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Figure 1: Study area 


Climate 

The climate of the study area is humid tropical type having an annual rainfall of about 2500 mm (at the coast) to 1220 mm (at the 
northern limit) which runs from April to October and characterized by a wet and dry season (Adebekun, 1978). The minimum and 
maximum temperature is about 22.48°C to 31.24°C with an average yearly temperature of about 26.6°C. 


2. METHODOLOGY 


In this study, land suitability analysis was integrated with RS and GIS in MCDA to determine the suitability of the study soil for rice 
production. The present study employed the following steps in achieving the methodology: 

1. Identification and development of a GIS database using multi-criteria suitability requirements. 

2. Processing and creation of RS and GIS data to compile criteria maps for suitability analysis. 

3. Determination of suitability level and scoring scale based on suitability criteria evaluation. 

4. Application of AHP for the assessment of the relative importance of each criteria weights. 

5. Application of a weighted overlay method for estimation of the land suitability index (LSI). 

6. Evaluation and reclassification of land suitability criteria into land suitability classes. 

7. Identifying areas with most suitable for rice cultivation. 


Rainfall Analysis 

The present study utilized the annual rainfall and Standardized precipitation index (SPI) approach forrainfall variability assessment of 
the study area.A total of 10 years (2009 to 2018) rainfall data were collected at the Department of Agro-meteorology and Water 
Management, College of Environmental Resources Management, FUNAAB. The SPI and annual rainfall data collected were used to 
assess the wetland potential of the study area for rice production. Table (1) shows the SPI classification for any location-based using 
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precipitation records (Doesken, and Kleist 1993). The SPI described by Akinsanola and Ogunjobi (2014) and Adegoke and Sojobi 
(2015) was used in this study and expressed as: 


SPI = _— atasatent Equation 1 


Where X: The rainfall in each particular year, X: The mean rainfall in each particular year,o: The standard deviation of rainfall in each 
particular year. 
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Figure 2: Soil classification of the study area 


Table 1: Annual standardized precipitation index given by McKee, Doesken, and Kleist (1993). 


Classification SPI 
Near Normal -0.99 to 0.99 
Moderately wet years 1.0 - 1.49 
Moderately dry years -1.0 to -1.49 
Very wet 1.5 to 1.99 
Severely dry years -1.5 to -1.99 
Wet extreme > +2.0 
Dry extreme <- 2.0 
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Processing and creation of spatial data layers 


Modified Normalized Difference Water Index 

Remote sensing data provide spectral indices for extracting information from Earth's surface such as waterbodies, soil, and surface 
temperature, etc. (Yang, et al., 2003). The present study employed the use of a Modified normalized difference water index (MNDWI) 
to effectively assess the presence of waterbodies in the study area. The MNDWI was used to remove or suppress the vegetation and 
soil noise effects from the study area (Xu, 2005). To perform the MNDWI, Landsat 8 Operational land images (OLI) and Thermal 
infrared sensor (TIRS) of the year 2015 and 2020 was acquired and downloaded from the United States Geological Survey (USGS) 
website repository. The acquired Landsat 8 OLI was further subjected to simplified data processing and user image interpretation in 
ArcGIS 10 ESRI (Environmental Systems Research Institute). Thereafter; a raster math calculator in spatial analyst tools was used to 
perform the MNDWI operation in the study area using the equation described by Xu, (2005): 


wi OITS (Green — MIR) . ae 
= (Green FR) quation 


Where: Green represents band 3, and MIR (Middle Infrared) represents and 6 of the Landsat 8OLI dataset. 


Enhanced Built-up and Bareness Index 

Application of Enhanced built-up and bareness index (EBBI) depends on the types of the satellite image, accuracy level, surface 
characteristics, and the purpose of the study (As-syakur et al. 2012). Study conducted by Li et al. (2017) revealed that EBBI serves as 
an efficient approach used in extracting built-up area and bareland cover information. Therefore, in this study, the EBBI was used to 
extracts information on the built-up and bare land using Landsat 8 dataset of the years 2015 and 2020. The EBBI of the study area 
was derived using the equation described by As-syakur et al. (2012): 


pper= (SWIR — NIR) eee 
(10 «root SWIR+TIR) quation 


Where; SWIR represents Short-wavelength infrared of band 6 in Landsat 8 OLI, 
NIR represents the Near-infrared of band 5in Landsat 8 OLI and 
TIR represents Thermal infrared of band 10 in Landsat 8 OLI. 


Vegetation cover analysis 

Vegetation indices are one of the primary sources of information used in monitoring the Earth’s vegetation cover (Gilabert et al. 
2002). The present study, utilized the Soil adjusted vegetation index (SAVI) to assess the vegetation cover of the study area using 
Landsat 8 datasets of the years 2015 and 2020. The SAVI was used due to its ability to correct the soil brightness in a study 
especially when the vegetation cover is low (Jiang et al. 2006). The SAVI described by Huete, (1988) and Verbesselt et al. (2006b) was 
used in this study and derived using the following equation: 


SAVI = eee) Equation 3 
Wteiepdap quation 


Where: NIR is the Near-infrared regions of band 5, Red is the visible red of band 4 and L is the constant or correction factor, ranges 
from 0 to 1. 


Land surface temperature 

Accurate Land surface temperature (LST) retrieval from remote sensing data depends on emissivity and geometry (Gao et al. 2013). 
For this study, LST was calculated using temporal information from the Landsat 8 OLI satellite images of the years2015 and 2020. 
The LST of the study area was subjected to two procedures. 


Step 1: The SAVI data derived from equation (3) were used to calculate the Proportion of vegetation (PV)of the study area which can 
be expressed as follows: 
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_ SAVI — SAVI min y 
~ ASAVI max — SAVImin 


vee ee es EQuation 4 


Where: SAV min was retrieved from the SAVI minimum value, which is -1, and SAVImax, retrieved from the SAV] maximum value, which 


is 1. Therefore, the PV of the study area was further expressed as follows: PV = ee) onnateaBQUuation.S 


Step 2: The thermal bands were converted to Digital numbers (DNs) to estimate the spectral radiance of the study area using bands 
10 and 11 in the Landsat 8 OLI image (Rasul et al.2015). The spectral radiance of the bands was described below and subjected to 
equation (6): 


LA = (0.0003342 * Band10 + 0.1 and LA = 0.0003342 * Band11 + 0.1)...... Equation 6 
LA = ML+QCAL + AL 


Where: LA is the TOA (Top of atmosphere) spectral radiance at the sensors aperture 
ML is the band-specific multiplicative rescaling factor from the metadata, 
Qcacis the quantized and calibrated standard product of the digital pixel values, 
AL is the band-specific additive rescaling factor from the metadata. 


Thereafter, Land surface emissivity (LSE) and Brightness temperature (BT) of the study area was calculated and described using 
equations (7) and (8). Hence, the LST of the study area was further calculated using equation (9) and expressed according to Jesus 
and Santana (2017): 


LSE = 0.004 * PV + 0.986... ................ Equation 7 
K2 
BT = —aq) _~ 273.15...» Equation 8 


Where: BT is the satellite brightness temperature (Celsius) 
K2 is the calibration constant 2 (Kelvin), thermal conversion constant from the metadata; 


K1 is the calibration constant 1 (Kelvin), thermal conversion constant from the metadata; 


BT 


LST = RRO Equation 9 


1 (=? *In LSE) 


Principal component analysis 

Principal components analysis (PCA) is a set of multivariate statistics that correlate variables and transform them into uncorrelated 
variables in the satellite image bands (Weng, 2012). The present study employed the use of PCA to improve the quality of the land 
use/land cover (LULC) classification of the study area (Singh and Harrison, 1985). To perform the PCA of the study area, Landsat 8 
OLI satellite of the year 2015 and 2020images were used in the ENVI software environment. The Landsat 8 satellite images were 
further subjected to data processing such as data redundancy, pre-processing, geometric adjustment, and conversion of DNs to 
reflectance values using Dark object subtraction (DOS) in ENVI software °'5 (Chavez, 1998; Chuvieco, 2002; Weng, 2012). The 
essence of this operation is to enhance the visual interpretation of the Landsat images for better enhancement during the analysis of 
the LULC change of the study area. The characteristics of the downloaded Landsat satellite images used for this study were shown in 
Table (2). The PCA of the study area was calculated and expressed according to Estornell et al. (2013) using equation (10).To 
calculate the eigenvector from the vector-matrix, an equation derived from Estornell et al. (2013) in equation (11) was used: 


1 1 

2 wit wie "2 
ye=[4 JSh 3 % : (| csv es esses ae ene Equation 10 

ee Wb,! Wb, b 2 


Where: Y is the vector of the principal components, w: is the coefficient of the transformation matrix of the eigenvectors in the 
diagonal covariance matrix of the original bands, and x: is the vector of the original data. 
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(Co = Ay Wy = 0 wee eee cee ee ee ee Equation 11 


Where: C is the covariance matrix,a,: is the k eigenvalues (Eleven in a Landsat OLI image composition), |: is the diagonal identity 
matrix, and Wx: is the K eigenvectors. 


Table 2: Characteristics of Landsat 8 OLI images used in the study 


Metadata 2015 2020 
Acquisition date 2015-10 - 12 2020 - 10 - 20 
Scene center time 10:03:03 10:03:05 

Path 191 191 

Row 055 055 

Sensor and platform OLI-TIRS OLI-TIRS 
Processing level Calibration Calibration 
Number of bands 11 11 

Sun azimuth 142.77069460 134.88500186 
Sun Elevation 52.10216187 50.77679746 


Land use/Land cover Analysis 

The acquired satellite images were enhanced in the Idrisi selva Software environment via (3 by 3) majority filter techniques for better 
visibility during analysis. True colour composite (TCC) was generated using suitable combinations of bands in the PCA images 
(d‘Entremont and Thomason, 1987; Good and Giordano, 2019). Considering the "Nigeria land use classification System" and the goal 
of this study, Anderson et al. (1976) land use classification scheme II and prior knowledge of the study area for over 6years was used 
to identify the Area of interest (AOI) features of the study area. For this study, an unsupervised and supervised classification 
algorithm was performed in the Idrisi selva environment to assess the nature of change in the study area using Landsat 8 OLI of the 
year 2015 and 2020 images. Among the supervisedclassification methods, the maximum likelihood algorithm is one of the most 
applied methodologies used in monitoring LULC changes especially in the field of agriculture (Biro et al. 2013; Zhang et al. 2015). 
Besides, maximum likelihood algorithm was considered due to its final outcome result and geovisualization (Silva and Costa, 2015: 
Liu, 2005; Sun et al. 2013). Therefore, this study employed the maximum likelihood classification system (MLCS) and the 
unsupervised K- means algorithm to assess the LULC change of the study area. The K- means of unsupervised classification were 
supported with MLCS due to the heterogeneous nature of wetlands (Harvey and Hill, 2001). The present study classified the LULC 
into 5 classes based on the MLSC algorithm techniques using the Landsat 8 OLI of the year 2020 (Table 3). The K- means algorithms 
of unsupervised classification were used to assess the Landsat 8 image of the year 2015. However, due to the complex nature of 
wetlands, such as seasonal variations and distinct landcover types, Google Earth's image was downloaded from the Google Earth 
engine repository website to validate the identified AOI of the study area. Thereafter, the classified LULC images of the study area 
were further subjected to post-classification processing. The essence of this operation was used to remove the misclassified pixels in 
the outlying LULC images in the study area. Shuttle radar topographic mapper (SRTM) of a 30-meter resolution was used to 
produced the Digital elevation model (DEM), aspect, stream order, and slope map of the study. The river and road map was also 
produced and used in the efficient mapping of the study area using ArcGIS 10 ESRI (Fig. 3a, 3b, and 3c). 


Table 3: Description of land use and land cover categories 


No Class name Description 

1 Built-up Lecture hall, offices, and laboratories 

2 Vegetation Mixed forest and grass 

3 Bareground Vacant land, open space, sand, bare soils, and landfill site 
4 Farmland Fallow land rainfed cropping planted cropping areas 

5 Wetlands Waterbodies, river, lakes, ponds creeks, and streams 
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Figure 3a: Slope and Aspect of the study area 
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Figure 3b: Slope and Aspect of the study area 
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Figure 3c: Slope and Aspect of the study area 


Precision Accuracy Assessment 

Precision accuracy assessment of the study area was carried out by comparing the extent of the LULC classes in the classified image 
to the reference data-set using the classified map (Congalton, 1991). The accuracy of the LULC map of the study area was computed 
through the following error matrix: Overall accuracy, kappa coefficient, producer's and user's accuracies (Foody, 2002). The overall 
accuracy is given by the ratio of the proportion of the correctly classified pixels to the total number of pixels in the error matrix © 
(Pontius Jr and Millones, 2011). To compute the error matrix in this study, the following formula was used. 


Area properly identified in a classification method 


Producer accuracy = ... Equation 12 


The area in the reference ground truth 


Area properly identified in a classification method 


User accuracy = ———— OOo Equation 13 


Total area calculated from the method 


Observed accuracy — Expected agreement 


Kappa coef ficient. = AADDMNMMWJ4____—+——_.. ... ... Equation 14 


1 — Expected agreement 


Soil Mapping and Analysis 

Reconnaissance and stratified random sampling were employed based on the soil survey manual (FAO, 2007a). A hard copy of the 
study area soil map was downloaded in the Geographic tagged image file (GeoTiff), scanned, and georeferenced in ArcGIS 10 ESRI 
(See figure 2). Coordinates of the study area soil types were retrieved from the scanned georeferenced soil map and further stored 
using a handheld Global positioning system (GPS) device. This operation was performed to allow accurate mapping and 
identification of the soil types in the study area. Thereafter, 5 representative soil samples were collected using a soil auger from each 
soil type at a predetermined interval making a total of 35 soil samples. The representative's soil samples were collected at the soil 
depth of 0 — 30cm for soil analysis. The soil samples collected were processed in the laboratory after air-drying at room temperature 
for soil analysis. The entire soil auger observations were stored and inputted with the help of handheld GPS for spatial mapping of 
the study area. Other physiographic features coordinate points were taken with the aid of a GPS and properly documented to assist 
in the adequate mapping of the study area. 
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Laboratory Data Analysis 
The following soil parameters were determined: nitrogen, potassium, phosphorus, hydrogen ions, soil texture, and soil pH. The 


particle size was determined by hydrometer method Bouyoucous described in Beretta et al., (2014). Total nitrogen was determined 
calorimetrically. Available phosphorus was extracted by Bray 1 method (Anderson and Ingram, 1993; Flavio et al., 2011) and 
phosphorus concentration was determined using a ultra violet spectrophotometer. Soil organic carbon was determined by 
dichromate oxidation procedure and organic matter was determined by 1.724 of carbon factor. Exchangeable bases (Calcium, 
Magnesium, Potassium, and Sodium) were extracted with neutral ammonium acetate. Ca and Mg were determined with Atomic 
Absorption spectrometry while K and Na were determined by flame photometer. Base saturation and effective cation exchange 
capacity were calculated. Soil pH was determined in a 1:2 soil to water suspension using a glass electrode pH meter. Heavy metal 
analyses, sub-samples (0.5 g) of each of the soil were digested. Digestion was done with 10 ml of a mixture of nitric (HNO3) and 
perchloric (HCIO4) acid in ratio 2:1 (v/v) for 90 min, initially at 150° C. After which 2ml of concentrated Hydrochloric Acid (HCL) was 
added to the mixture. The temperature of the digest was then increased to 230° C for another 30 minutes on the digester. On 
completion of digestion, digests were allowed to cool down at room temperature. Thereafter, the content of each digestion tube 
was transferred into a 50 ml volumetric flask and made to volume with distilled. 


Data analysis 

Descriptive statistics and geostatistics were used to analyze and interpret the soil datasets and spectral index. Descriptive statistics 
were run in the R package and the pairwise Pearson's correlation between the spectral indices was calculated in RStudio v4.0.3. 
Thereafter, a correlation plot was generated using the R package “corrplot.” Coordinates of the analyzed soil data were saved in 
excel spreadsheets as text delimited and each of these soil data concentrations was plotted and display in the ArcGIS 10ESRI. 
Thespatial distribution of collected soil nutrients was prepared using Inverse distance weighted (IDW) techniques for soil data 
interpolation in ArcMap (Li and Heap 2008). IDW is a technique that relies on the significance of known points from the output point 
(Watson and Philip 1985). Therefore, in this study, the equation (15) derived by Agris (1998) was used to interpolate the estimated 
values and the weight of the collected soils sample was established using equation (16): 


Z * (xo) = Vp wiZ (xi) ...........Equation 15. 


Where: wi, the weight assigned to the value at each location: Z (xi) and n, number of close neighboring sample data points used for 
estimation. Thereafter, a skewness value of < 1 was employed using the IDW (Agris 1998). 


; 1/dpi 
wee /dp 


i=1 


Where: dp denotes the distance between point i and the unknown point; p. exponent parameter. 


Determining criteria weights using AHP 

AHP is a method that combines two or more criteria at a time through a pairwise comparison matrix and consistency ratio (Saaty, 
1987). The AHP relies on the relative importance of weight given to one criterion over another (Saaty 1980). The present study 
utilized the spectral indices (LST, MNDWI, EBBI, SAVI), classified LULC, and the interpolated soil data to determine the relative 
importance weight of each criterion. The relative preference of the criteria was subjected to a pairwise comparison scale (Table 4). 
MCDM was used to produce the criteria maps using the calculated weights for each criterion as attribute data. The pairwise 
comparison matrix approach and consistency ratio (CR) was computed as: 


(Amax= 1) 


CR= 7] bes ses sav eet ees tas sss ees ses ses asses ase BQuation 17 


Where Amax : is the principal eigen value, n: is the number of elements compared, and RI: is the random consistency index value that 
depends on the number of criteria that are being compared. The overlay weighted method in spatial analyst tool of ArcGIS 10ESRI 
was used to estimate the LSI of the study area according to the formula described by Saaty, (1987): 
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Where Vi: is the suitability index in each raster cell, w;: is the weight of criterion, j and n is the total number of criteria. The raster layer 
of the resulted map was further classified into the adopted FAO (1976) following the procedure described by Sys et al. (1991) land 
suitability evaluation (Table 5). Thereafter, the classified land suitability map was overlaid with the soil map of the study area. 


Table 4: Scale for pairwise comparisons (Saaty 1980) 


Intensity of Importance 


Definition 


OON DWN KRWDY | 


Equal importance 


Equal to moderate importance 


Moderate importance 


Moderate to strong importance 


Strong importance 


Strong to very strong importance 


Very strong importance 


Very to the extreme importance 


Extreme importance 


Table 5: Land suitability class description and the corresponding suitability index 


Class Suitability 


Suitability Index 


Description 


S1 Optimally suitable 

S2 Moderately suitable 
S3 Marginally suitable 

N1 Marginally not suitable 
N2 Permanent Unsuitable 


> 75 
50 - 75 


25 - 50 


10-25 


< 10 


Land without any significant limitation 

A moderately severe limitation which 
reduces productivity 

Overall severe limitations; given land 

use is only marginally justifiable 

The land-use type under analysis is not 
acceptable at all for the land 

The land-use type under analysis is 
permanently not acceptable at all for the 
land 


3. RESULTS AND DISCUSSION 
Rainfall Analysis 


Table 6: Classification of annual rainfall based on SPI 


Classification 


Years 


Dry extreme Years 
Moderately dry years 
Moderately wet years 
Near normal years 
Very wet years 

Wet extreme years 


2009, 2013, 2017, 2015, 2018. 


2010, 2011, 2012, 2014, 2016. 


Rainfall data analysis and SPI calculation were used to assess the dry and wet years of the study area. Classification of the annual 
rainfall was based on the SPI of the study area. The result showed that the SPI of the studied area falls into moderately wet and very 


wet years (Table 6). The present result suggests that the precipitation of the studied area can supports wetland potential for 


optimum rice production. 


Analyses of LST, MNDWI, EBBI and SAVI of the study area 
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Assessment of wetlands potential using spectral satellite image data is widely used due to its ability to provides higher accuracy and 


tends to be less time-consuming ’® lensen, 1996). In this study, spatial variation of extracted spectral indices of LST, MNDWI, EBBI, 
and SAVI has been studied for a period of 5 years from 2015 to 2020 (Figures, 4 to 7). The spatial variation of the LST for both years 
elaborates that in the year 2015, minimum and maximum temperature experienced in the studied area ranged from 23.4° K to 29.5° 
K (Celsius). For the year 2020, the maximum temperature tends to decrease from 29.5° Kto 25.1° K and the minimum temperature 
remains 23.4° K in the studied area (Figure 4).According to Sys et al. (1993), LST that ranges from 23.4° K to 25.1° K can be 
considered favorable for optimum rice production. Also, Liu et al. (2016) and Wang et al. (2018) posits that area or region with low 
or higher LST is one of the crucial environmental problems faced in the field of agriculture especially wetland soils. As seen in 
Figures5, 6, and 7,it can therefore be inferred from the resulted images that MNDWI depicts higher waterbodies coverage in the year 
2020 than 2015. McFeeters, (1996) and Chowdary et al. (2010) noted that MNDWI serves as one of the well-established means of 
detecting surface waterbodies due to higher reflectance in the SWIR compared with the NIR wavelength range. Similarly, the EBBI in 
the present study was quite high in the year 2020 indicating the percentage of the built-up and bare land area increased more than 
EBBI in the year 2015. Therefore, the EBBI used in this study reveals a better accuracy to the built-up areas with a higher degree of 
correlation in assessing LULC change. The research conducted by Sinha et al. (2016) and Bouhennache et al. (2015) also 
corroborated with the present study. Besides, the increased experience in both EBBI and MNDWI in the year 2020 could also be 
traced to anthropogenic activities such as reduction or loss of vegetation cover, etc. Considering the result of the SAVI, it shows that 
the study area demonstrates very low vegetation densities in the year 2020 than in 2015. The decrease in the year 2020 SAVI could 
also be responsible for an increase in the EBBI in the year 2020 of the study area. Moreover, Aman et al. (1992) opined that loss or 
low vegetation directly increases human activities, especially in developing countries such as Nigeria. The MNDWI, EBBI, LST, and 
SAVI used in this study serves as indispensable spectral indices to assess the rate of change in the study area (Yang, et al., 2003). 
Tables (7) shows the descriptive statistics of the spectral index of the study area. The statistical data result of the LST, MNDWI, EBBI, 
and SAVI was assessed and recorded from each indices equation. 


Table 7: The statistics summary of each index used in the study area 


Statistic MNDWI EBBI SAVI LST 

2020 2015 2020 2015 2020 2015 2020 2015 
Mean - 0.05 - 0.20 - 0.31 - 0.22 0.45 0.36 26.40 29.66 
Std.dev 0.48 0.20 0.88 0.82 0.25 0.15 1.76 2.97 
Min. - 0.61 - 0.72 - 1.88 - 1.38 0.03 0.05 23.43 25.13 
Max 1.20 0.09 1.32 1.25 1.01 0.57 29.50 35.42 


Min: Minimum; Max: Maximum; Std.dev: Standard deviation. 
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Figure 5: Modified Normalized Difference Water Index of the study area 
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Figure 6: Enhanced Builtup and Bareness Index of the study area 
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Figure 7: Soil Adjusted Vegetation Index of the study area 


Correlation between spectral indices 

The spectral index of LST with MNDWI (r? = 0.80, P< 0.02), and LST with SAVI (r* = 0.73, P< 0.04), had a strong correlation while LST 
with EBBI (r? = 0.58, P< 0.08) had a moderate correlation. In general, a strong association between the spectral indices was observed 
(Figure 8). 
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Figure 8: Correlation matrix of the study area. 


Principal Component Analysis 

The present study employed the principal component techniques to improve the quality of the LULC classification of the study area. 
In this study, the PCA was applied to the seven bands in the satellite images resulting in the eigenvalues; mean; standard deviation; 
minimum, and maximum (Tables 8 and 9). The seven bands of satellite images used in this study were employed to remove data 
redundancy and also to save computation space and time during analysis. As seen in Tables 10 and 11, the PC1, PC2, and PC3 in the 
year 2015were negatively loaded on the bands 5 and 6 in the following order: -0.481, -0.476, -0.381, -0378, and -0.257. It was also 
observed that PC3, PC4, PC6, and PC7 shows a high positive loading on band 1 (0.623), band 5 (0.683), band 6 (0.813), and band 7 
(0.680). In contrast, the PCA of the year 2020 had a negative loading in almost all the image bands. For the positive loading,PC4, 
PC5, and PC6 show a higher positive loading in band 3 (0.703), band 2 (0.963), and band 4 (0.560). The result of the PCA in this study 
denotes that satellite bands that are positively loaded mainly represent no change in the overall brightness. The negatively loaded 
measure the changes in brightness of the two images used in the present study (Koutsias et al. 2009). Also, the application of PCA in 
this study serves as one of the important pre-processing techniques used in interpreting LULC change in a study (Rojas 2011). 


Table 8: The Year 2015 PCA of the study area 


Bands Minimum Maximum Mean Standard Deviation Eigen Values 
PCI 0.013780 0.020317 0.014558 0.000478 0.000025 
PC 2 0.011631 0.019592 0.012552 0.000639 0.000009 
PC 3 0.009226 0.020314 0.010804 0.000875 0.000001 
PC 4 0.006932 0.024129 0.009045 0.001593 0.000000 
PC5 0.008665 0.042075 0.027242 0.002975 0.000000 
PC6 0.004266 0.046896 0.019923 0.003177 0.000000 
PC7 0.002324 0.051162 0.010936 0.003393 0.000000 


PC: Principal component. 


Table 9: The year 2020 PCA of the study area 


Bands Minimum Maximum Mean Standard Deviation Eigen Values 
PC 1 0.018748 0.020297 0.019207 0.000210 0.000009 
PC 2 0.017424 0.019254 0.017896 0.000271 0.000001 
PC 3 0.015415 0.018353 0.016065 0.000411 0.000000 
PC 4 0.014264 0.019298 0.015470 0.000721 0.000000 
PC5 0.018575 0.028806 0.024057 0.000836 0.000000 
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PC6 0.013239 0.053289 0.022143 0.001866 0.000000 
PC7 0.007797 0.135582 0.013936 0.002293 0.000000 
PC: Principal component. 
Table 10: The Year 2015 PCA Eigenvector of the study area 
Eigenvector PC 1 PC 2 PC 3 PC 4 PC 5 PC6 PC7 
Band 1 0.084 0.116 0.162 0.309 -0.081 0.623 0.680 
Band 2 0.013 0.019 -0.016 0.043 -0.984 -0.172 0.019 
Band 3 0.262 0.320 0.405 0.479 0.136 -0.620 0.182 
Band 4 0.157 0.198 0.287 0.394 -0.071 0.436 -0.709 
Band 5 -0.481 -0.476 -0.257 0.683 0.034 -0.089 -0.022 
Band 6 -0.378 -0.381 0.813 -0.224 -0.034 -0.001 0.017 
Band 7 0.725 -0.688 0.025 0.041 -0.002 0.000 0.002 
PC: Principal component. 
Table 11: The year 2020 PCA Eigenvector of the study area 
Eigenvector PC 1 PC 2 PC 3 PC 4 PC5 PC 6 PC 7 
Band 1 -0.055 -0.075 -0.116 -0.218 -0.058 -0.603 -0.751 
Band 2 -0.006 0.003 0.047 0.006 0.963 0.159 -0.211 
Band 3 -0.156 -0.214 -0.333 -0.420 -0.173 0.703 -0.344 
Band 4 0.243 0.283 0.354 0.560 -0.190 0.339 -0.521 
Band 5 -0.725 -0.390 -0.109 0.556 -0.007 -0.019 -0.037 
Band 6 -0.596 0.429 0.570 -0.363 -0.039 0.037 -0.009 
Band 7 0.183 -0.729 0.641 -0.149 -0.033 0.020 -0.011 


PC: Principal component. 


Land-use changes from 2015 to 2020 

The analyzed TCC (True color composite) reveals that there are 5 major types of LULC in the study area namely: bareground, 
farmland, wetland, vegetation, and built-up areas. The land-use classification maps obtained from the TCC was displayed in Figure 9. 
The result of the land use classification maps showed that wetland and farmland are the dominant land use classes with total areas 
of about 60% (Table 12). It was observed from the land use classification map that built-up, wetland, bareground, vegetation and 
farmland had respective areas of 2.7%, 22.6%, 19.7%, 22.8%, and 32.3% in the year 2015. Consequently, in the year 2020, vegetation 
decreased from 2277.1 ha to 1098.3ha, and an increase was observed in farmland, wetland, bareground, and built-up areas in the 
following order: 33.9%, 27.8%, 24.0%, and 3.3% respectively. The increase observed in the wetland in the year 2020 could contribute 
to regulations of climate change in the study area. Mitsch et al. (2012) opined that an appreciable increase in wetlands serves as a 
climate change regulator and supports atmospheric carbon-dioxide reduction. The result of the post-classification showed that the 
area under wetland, bareground, farmland and built-up areas increased in the following order: 517.9ha, 431.8ha, 169.2ha, and 
59.9ha, and a decrease was observed in vegetation cover (Table 13). The decreased in vegetation observed may be traced or 


attributed to an increase in human activities in the study area during the past 5 years. 
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Figure 9: Land use and land cover of the study area 


Table 12: Land use and land cover classification of the study area. 


LULC class 2015 2020 
Ha % Ha % 

Builtup 266.7 27, 326.6 3.3 
Wetland 2263.3 22.6 2781.2 27.8 
Bareground 1972.9 19.7 2404.7 24.0 
Vegetation 2277.1 22.8 1098.3 11.0 
Farmland 3228.2 32.3 3397.4 33.9 
Total 10008.2 100 10008.2 100 


Ha: Hectare, %: Percentage 


Table 13: Change detection of the study area 


LULC class Area (Ha) in year Area (Ha) in year Change detection (Ha) 
2015 2020 between 2015 and 2020 

Builtup 266.7 326.6 59.9 (Increase) 

Wetland 2263.3 2781.2 517.9 (Increase) 

Bareground 1972.9 2404.7 431.8 (Increase) 

Vegetation 2277.1 1098.3 1178.8 (Decrease) 

Farmland 3228.2 3397.4 169.2 (Increase) 


Ha: Hectare 


Accuracy Assessment and Analysis 

Accuracy assessment is the most important part of remote sensing and land-use change classifications (Li, 2005). The kappa 
coefficient and overall accuracy were used to verify the accuracy of the land use classification of the study area. In this study, the 
classification overall accuracy efficiency reached 89.47% and 92.22% respectively. The MLSC algorithm applied in the Landsat 8 
satellite images in this study showed that wetland, farmland, vegetation, background, and builtup areas were accurately classified 
with 83.33% and 88.57% kappa coefficient (Tables 14 and 15). 


Table 14: Error matrix and accuracy assessment of the land use map for 2015 
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Producer User 


Land use Wetland — Farmland Builtup Bareground Vegetation Total Accuracy Accuracy 

(%) (%) 
Wetland 40 0 0 10 3 53 93.02 75.47 
Farmland 0 72. 0 0 3 75 93.51 96.00 
Builtup 1 3 31 1 3 39 86.11 79.49 
Bareground 2 2 0 60 0 64 84.51 93.75 
Vegetation 0 0 5 0 52 57 86.25 91.23 
Total 43 77 36 71 61 288 


Overall accuracy = 89.47 % Kappa coefficient=83.33% 


Table 15: Error matrix and accuracy assessment of the land use map for 2020 


Producer User 


Land use Wetland = Farmland Builtup Bareground Vegetation Total Accuracy Accuracy 

(%) (%) 
Wetland 30 0 2 0 10 42 78.94 71.42 
Farmland 7 59 0 5 3 74 85.50 79.72. 
Builtup 0 0 92 12 0 104 97.87 88.46 
Bareground 0 10 0 112 0 122 86.82 91.80 
Vegetation 1 0 0 0) 300 301 95.84 99.67 
Total 38 69 94 129 313 643 


Overall accuracy = 92.22 % Kappa coefficient = 88.57 % 


Descriptive Statistics of the analyzed soil properties 

The descriptive statistical summary of the analyzed soil properties was presented in Table (16). The soil pH values were ranging from 
5.2 to 8.63 with a mean of 6.71, which was also similar to the median values of 6.61. The concentration of soil Organic matter (OM) 
ranged from low to high (0.6 to 5.8%), with a mean of 3.22%. Total nitrogen (TN) was relatively low (ranging from 0.29 to 0.03%), 
with a median of 0.17%, and a mean of 0.16%. Available phosphorus, potassium, and Cation exchanged capacity (CEC) were within 
their respective medium ranges. The micronutrients measured showed that zinc was very high (range > 1.0 mg kg"!) with a median 
of 9.4 mg kg*'and a mean of 10.61mg kg"'whilecopper and magnesium were also very high. For copper, the range was> 0.5mg kg” 
with a mean of 2.98mg kg™!anda median of 2.40mg kg~'. Magnesium also had a mean of 1.64mg kg™'with a median of 1.57mg kg” 
‘It was observed that the majority of the soil properties had low skewness with a soil texture ranging from 4.58 to 90.5 % for sand, - 
0.39 to 12.5% for silt, and 1.88 to 17.0% for clay. °*Kravchenko and Bullock (1999) reported that soil properties with low skewness 
(<1) yielded the most accurate estimates of soil properties. The present study also interdem with the work done by Pal et al. (2010); 
Sharma et al. (2011); Wang and Shao (2013) that estimating spatial variability of physical and chemical soil properties serves as one 
of the prerequisites for soil and crop-specific management techniques. 


Table 16: Descriptive statistics for the studied soil properties 


‘ ; . ; Standard 
Variable Min Max Mean Median Skewness Kurtosis as 
Deviation 
Soil pH 5.20 8.63 6.71 6.61 0.56 1.14 0.71 
Organic matter 0.69 5.85 3.22 3.34 0.09 29 1.70 
Total Nitrogen 0.29 0.03 0.16 0.17 0.05 -1.26 0.08 
Phosphorus 0.15 18.52 2.47 1.41 3.27 10.54 4.18 
Potassium 0.06 0.36 0.18 0.18 0.49 0.70 0.06 
Zinc 7.00 25.00 10.61 9.4 2.77 7.51 4.25 
Magnesium 0.72 2.71 1.64 1:57. 0.29 0.05 0.47 
Copper 0.90 10.4 2.98 2.40 2.83 9.86 1.88 
Sand 74.5 90.5 84.50 86.5 -0.74 -0.44 4.58 
Silt 0.5 12.5 6.74 6.5 0.20 -0.39 3.21 
Clay 7.00 17.0 8.76 8.00 1.88 4.18 2.45 
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CEC 3.99 17.37 10.32 9.47 0.61 0.14 3.45 
Min: Minimum; Max: Maximum; CEC: Cation exchanged capacity. 


Digital soil mapping using IDW 

The variability of soil nutrients across the study area calls for a proper soil management strategy. The spatial distribution map of the 
selected chemical soil properties was produced by IDW (Inverse distance weighted) techniques and the results were grouped into 
various classes based on the range representing their magnitude in the studied soil (Figures 10 to 12). 
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Figure 10: Soil pH and Soil OM 


Soil pH varied from strongly acidic (< 5.5) to moderately alkaline (>8.4) in the study area. These results are in agreement with 
those reported in a recent study by Tobore et al., (2020). However, the variation in the soil pH observed could be attributed to the 
nature of the alluvial parent material, micro-topography, and the types of fertilizer used at the time of cultivating the study area °” 
(Vasu et al. 2017). According to Bremner and Mulvaney (1982), Westarp et al. (2004), Regmi and Zoebisch (2004) reported that soil 
management practices such as crop nutrient uptake and harvest without replenishment with poor crop residues management could 
be responsible for the acidic nature in the present studied soil. 


www.discoveryjournals.org | OPEN ACCESS 


pase | 9 


ARTICLE 


3°20'0"E 3°23'20"E 3°26'40"E z 3°20'0"E 3°23'20°E 3°26'40"E z 
r-4 zo ra 
° OR a 
S aR R 
KR RK 
z < 
z z3 = 
8 2 © 
© OR KR 
KR KR 
z = 
z za g 
a ee 2 
2 2K R 
KR RK 
Legend 
Legend 
Study Boundary (FUNAAB) 
Study Boundary (FUNAAB) ® Soil Sampling 
© — Soil Sampling 2 || Total Nitrogen é 
z ‘ zo io] 
= || Available Phosphorus 5 27) (%) =) 
21]| (mg/kg) Sk < 0.1 (Very Low) f 
~]| (] <15 dow = HE 01-02 (Low) 
HE 15-25 (veaium) HEE 0 2-05 (medium) 
3°20'0"E 3°23'20"E 3°26'40"E 3°20'0"E 3°23'20"E 3°26'40"E 
Figure 11: Total N and Available Phosphorus 
3°20'0"E 3°23'20°E 3°26'40"E > 3°20'0"E 3°23'20"E 3°26'40"E z 
ca zo ° 
° OR a 
a Qik KR 
KR RK 
=e = 
Zz z¢ 2 
S 26 © 
© OR R 
BS KR 
= = 
R Re ® 
2 as . 
KR KR 
Legend 
9 Legend 
Siidy Soundary (FUNAAB) Study Boundary (FUNAAB) 
2 Sorisampiny @ = =Soil Sampling 
Cation Exchanged Capaci Zz Zz 
a 9) pacity = 5 || Potassium S 
£ cmol/kg © 21} cmol/kg = 
21} | <8 (ery Low) 2 R |] 7) <02 Wer tow) R 
~ || EB 6-12 00w ie HM 02-03 (Low) 
HE 12 - 25 (voderate) km || BB 0.3 - 0.6 (Moderate) 
3°20'0"E 3°23'20°E 3°26'40"E 3°20'0"E 3°23'20"E 3°26'40"E 


Figure 12: CEC and Potassium 


The soil OM (Organic matter) was relatively low (<2%) in the majority of the study area, followed by medium (2 — 3%) and high in 
> 3%. The low organic matter content in the study soils could be traced to the acidic nature of the soil pH thereby leading to low 
levels of soil OM in the studied soil. Besides, poor management practices and loss of vegetation could also be responsible for the 
low soil OM in the studied soil. 
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Total nitrogen (TN) was deficient in the study area with values < 0.2 % (low and very low) and medium with > 0.2 - 0.5%. The low 


variation in TN content in the study soil may be related to low soil OM content and the application of inappropriate fertilizer. For 
instance, the increased rate of mineralization and insufficient application of nitrogen fertilizer to crops such as maize could also be 
responsible for the low of TN in the studied soil (Sherchan and Gurung, 1995). 

The available phosphorus was medium with values > 15 and low with values <15 in the study area. The low soil OM may account 
for the low level of available phosphorus in the studied soils. 

The potassium concentration in the study soil ranged from very low (<0.2cmol/kg), low (0.2 — 0.3cmol/kg) and moderate (0.3 — 
0.6 cmol/kg) in the study area. The majority of the CEC was low with values ranging from very low (<6 cmol/kg), low (6 — 12cmol/kg), 
and medium (12 -25cmol/kg) in the study area. The effect of the low potassium and CEC may be attributed to continuous or yearly 
cultivation of the study soils for cropping practices. 


Determining criteria weights using AHP 

The present study employed AHP as a systematic approach to conduct MCDM in the study area. To produce the land suitability map 
of the study area for rice production, the criteria maps produced were used in line with the result of the AHP in MCDM. The criteria 
spatial layers of the relative importance were derived and the most important criterion for decision making used in this study was 
listed in Table (17) based on the pairwise comparison. The common scale of 0 to 100% was obtained from the AHP procedure and 
the consistency ratio of the pair-wise matrix was 7.5%. This implies that the comparison is less than 10% and thus considered 
suitable to weighted overlay using each criterion. 


Table 17: Resulting weights for the criteria based on pairwise comparison 


. Priority 
LST MNDWI LULC SAVI EBBI Soil data (%) ank 
LST 1 2 0.5 2 2 2 23 1 
MMDWI 0.5 1 1 1 0.5 21 2 
LULC 2 1 1 1 2 1 20 3 
SAVI 0.5 1 1 1 0.5 2 14 4 
EBBI 0.5 1 0.5 2 1 1 12 5 
Soil data 0.5 1 1 0.5 1 1 10 6 


MCDM within GIS 

The potential wetland suitability of the study area was generated by integrating the thematic map layers of LST, MNDWI, LULC, SAVI, 
EBBI, and analyzed soil data in ArcGIS 10 ESRI. The integrated raster layers were classified into three classes namely: marginally not 
suitable (NI), marginally suitable (S3), and moderately suitable (S2). The classified land suitability map was overlaid with the soil map 
of the study area. The result of the land suitability map shows that 28.9% of the area was moderately suitable, 33.2% was marginally 
suitable, and 37.9% was marginally not suitable for rice production (Figure 13). The result of the marginally not suitable (N1) could 
be traced to anthropogenic activities in the study area. The area with moderately (S2) and marginally (S3) suitable can be used for 
rice cultivation in the study area. However, there is a need to add appropriate fertilizers (from organic and inorganic sources) to the 
soil to increase fertility status. 


4. CONCLUSION 


The importance of wetlands to humankind and the environment requires spatial immediate intervention especially in developing 
countries such as Nigeria. The main objective of this study was to assess the wetland's potential suitability for rice production using 
remote sensing data and a GlS-based conceptual framework. The GIS-based conceptual framework combines through AHP has 
helped to adequately and efficiently map the study area for rice production in less time. The selection of multi-criteria factors used 
in this study shows that spectral indices (such as LST, LULC, MNDWI, EBBI, and SAVI) and the interpolated soil maps can help to 
facilitate easy identification of wetlands potentials and assist farmers in identifying the expected nutrient levels in the study area. 
Besides, the spatial distribution and mapping will aid farmers and various soil users to understand the present conditions of the 
different land use of the study area. These results can be implemented by the GIS techniques to overcome future threats that may 
pose to the wetland of the study area. 
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Figure 13: Land suitability for rice production. 
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